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ABSTRACT 

A variety of independent observational studies have now reported a significant 
decline in the fraction of Lyman-break galaxies which exhibit Lya emission over the 
redshift interval z = 6-7. In combination with the strong damping wing extending 
redward of Lya in the spectrum of the bright z = 7.085 quasar ULAS 1120+0641, this 
has strengthened suggestions that the hydrogen in the intergalactic medium (IGM) is 
still substantially neutral at z ~ 7. Current theoretical models imply H i fractions as 
large as 40-90 per cent may be required to explain these data assuming there is no 
intrinsic evolution in the Lya emitter population. We propose that such large neutral 
fractions are not necessary. Based on a hydrodynamical simulation which reproduces 
the absorption spectra of high-redshift (z ~ 6 — 7) quasars, we demonstrate that 
the opacity of the intervening IGM redward of rest-frame Lya can rise rapidly in 
average regions of the Universe simply because of the increasing incidence of optically 
thick absorption systems as the tail-end of reionisation is approached. Our simulations 
suggest that a volume-averaged neutral fraction of ~ 10 per cent is sufficient to explain 
the Lya emitter fraction at z = 7, assuming the observed decline is due to evolution in 
the IGM H i fraction only. These data do not require a large change in the IGM neutral 
fraction by several tens of per cent from z = 6-7, but may instead be indicative of the 
rapid decrease in the typical mean free path for ionising photons expected during the 
final stages of reionisation. 

Key words: dark ages, reionisation, first stars - galaxies: high-redshift - intergalactic 
medium - quasars: absorption lines. 



1 INTRODUCTION 



The opacity observed blueward of rest frame Lya in the 
spectra of distant quasars rises t oward higher redshifts 
jFan et al.ll2006l ; iBecker et ai]|2007h , indicating the fraction 
of neutral hydrogen in the intergalactic medium (IGM) is 
small but increasing with lookback time. This observation, 
coupled with the Thomson optic al depth measured from cos- 
mic microwave background data (|Komatsu et al.|[201ll k sug- 
gests that an extended epoch of hydrogen reionisation was 
ending by z ~ 6 — 7. Reaching further into the epoch of reion- 
isation - when intergalactic hydrogen is still substantially 
neutral - is difficult, but Lya selected galaxies may offer 
one promisi ng route to probing somewhat deeper into thi s 
distant era (|Miralda-Escude fc Reesi Il998l ; iHaimanl 120021} . 
The increasing Hi content in the IGM at z > 6 produces a 
Lya damping wing which can extend redward of a galaxy's 



Lya emission line. T he damping wing reduc es the visibility 
of the Lya emission (jMiralda-Escu dc 1998), and the num- 
ber of Lya emitting galaxies (LAEs) observed in flux-limited 
surveys will thus decrease as the ambient Hi fraction rises 
with increasing redshift. 

The resonant nature of the Lya transmissio n un- 
fortu nately complicates thi s simple picture (e.g. ISantos! 
12004) . Source clus tering (|Furlanetto et all 120061 ). dust 
|Hansen fc Ohl [20061 ). and resonant Lya scattering within 
the sources themselves - w hich is sensitive to both the 
H i content and gas velocities (|Diikstra fc Wvithj2010h - all 
impact on the visibility of Lya emission. The complex reso- 
nant radiative transfer within the galaxies often leads to sub- 
stantial line redshifts of several hundred kms -1 which can 
significantly enhance the visibility of LAEs when the sur- 
rounding IGM is still substantially neutral (|Diikstra et al.l 
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20 111 ; lLaursen et alj|201l|; IZheng et all l201ll ; iBarnes et all 
201 J ; IJeeson-Daniel et afll2012l) . 



Recently, however, a significant decrease in the trans- 
mission of Lya photons f rom high redshift galaxies has been 
found between 2 = 6 



20 111; lHaves et all 12011 



is from high rcdsmlt galaxies has beer 
7, llStark et al.ll2010iriPentericci et"al 



2012;; ICurtis-Lake et alj 



lOno et al, f |2012| ; I Schenker et al 
2012h . This result is based on the 



rapid decline in the fraction of Lyman break galaxies (LBGs 
which exhibit Lya emission - a quantity which is (in princi- 
pie) less susceptible than t he LAE luminosity function (e.g 



Malhotra fc Rhoadsl liooi ; iKashikawa et alj|2006l ; iHu et all 
20101 ) to observational selection effects and intrinsic evolu- 



tion in the LAE population. Based on comparisons to ex- 
isting numerical simulations and semi-numerical models of 
reionisation, the reported rapid decrease in the transmis- 
sion of Lya emission at 2 ~ 6 — 7 may be interpreted as 
a rapid change in the volume-averaged neutral fraction (by 
seve ral tens of percent) o v er a rather short redshift inter- 
val llPentericci et al J [201 it ISchenker et al.ll2012l ; lOno et all 
12012] ). 

Such a rapid change in the neutral fraction is, however, 
at odds with the rather low ionising e missivity suggested 
by the Lya forest data at z = 5-6 (|Bolton fc Haehneltl 
l2007bl ) which appears to imply a much slower evolu- 
tion of the neutral fraction. Theoretical reionisation mod- 
els matching the Lya forest da t a (e.g. Ciardi et al.l |2012| ; 
iKuhlen fc Faucher-Giguerd 120121 ; iJensen et al.ll2012T ) there- 
fore have considerable difficulty in explaining the rather 
large neutral fractions (~ 40-90 per cent) that have been 
inferred from the recent LAE/LBG observations at z ~ 7. 
Alternative explanations which - either individually or in 
combination - might allow for a more modest change in the 
IGM neutral fraction include an increase in the escape frac- 
tion of ionising photons, or (perhaps somewhat less plau- 
sibly) an increase in the interstellar d ust content of the 
galaxies toward higher red shift (see e.g. iDaval et alj l20lH 
iForero- Romero et alj|2012l) 

In this paper, we shall argue that one other clue to 
this puzzle is provided by the quasar ULAS Jl 120+0641 
at z = 7.085, recently discov ered by the UKIRT In frared 
Deep Sky Survey (UKIDSS, lLawrence et~all 120071 ). The 
spectrum of ULAS J1120+0641 exhibits a strong Lya damp- 
ing wing extending red w ard of the quasar Lya emission 
line |Mortlock et alj|201ll ) - exactly what is needed to sup- 
press the generally redshifted (relative to the systemic red- 
shift of the galaxy) emission of LAEs. Theoretical models 
predict that the regions surrounding rar e, bright quasars 
are a mongst the firs t to reionise (e.g . iFurlanctt o et alj 
20041; llliev et al.ll2006l; iTrac fc Cenl 120071 ; IZahn et alj|2007l ; 



McQuinn et all l2007bj ). If the environment around ULAS 
J1120+0641 is typical of quasar host galaxies at 2 > 7, 
the significantly weaker ionising radiation field expected 
around even the brightest LAEs implies the impact of the 
red Lya damping wing on LAE visibility should be even 

stro nger. 

iBolton et all |201ll ) used high resolution simulations 
of Lya absorption to demonstrate that unless a proxi- 
mate damped Lya absorber lies within ~ 5pMpc of the 
quasar, a volume averaged Hi fraction of at least (/hi)v = 
0.1 is required to reproduce the damping wing observed 
around ULAS J1120+0641. In this work we show that an 
H 1 fraction of not more than 10 per cent could also explain 



the rapid evolution of the LAEs - obviating the need for 
Hi fractions as large as 40-90 per cent. The difference in 
the inferred neutral fraction is due to the inclusion of small- 
scale (~ 20pkpc), optically thick absorptions systems which 
remain unresolved in large-scale ( > 100 cMpc) reionisation 
simulations. On including these systems the Lya opacity 
up to a few hundred kms -1 redward of rest frame Lya rises 
more rapidly with H 1 fraction than usually predicted. As we 
shall demonstrate, this is because the intervening Lya opac- 
ity of the IGM at the end of hydrogen reionisation is strongly 
affected by the increasing size of these systems, which nat- 
urally occurs as the metagalactic photo-ionisation rate de- 
creases. A more modest neutral fraction of < 10 per cent is 
much easier to reconcile with the fact that quasar absorption 
spectra at 2 ~ 6 appear to i ndicate that the IGM is highly 
ionise d only 180 Myr later dFan et all I2OO6I ; IWvithe et alj 
|200S| . but see also lMesingerll2010l ). 

We begin our analysis in Section 2, where we describe 
our modelling of the IGM Lya opacity and the procedure we 
adopt for incorporating self-shielded gas into a hydrodynam- 
ical simulation. In Section 3 the properties of the optically 
thick absorbers in the simulations are discussed in more de- 
tail. We explore the implications of these absorption systems 
for the visibility of Lya emission from galaxies in Section 4, 
and in Section 5 we conclude. We assume the cosmological 
parameters Cl m = 0.26, fi A = 0.74, Q b h 2 = 0.023, h = 0.72 
and a helium mass fraction of Y = 0.24 throughout. Comov- 
ing and proper distances are denoted by using the prefixes 
"c" and "p", respectively. 



2 NUMERICAL SIMULATIONS OF THE IGM 
2.1 Hydrodynamical model 

We model the IGM using a high-resolution hydrodynam- 
ical simul ation performe d with the parallel Tree-SPH code 
GADGET-3 (|Springel2005l ). The simulation employs a box size 
of Wh^ 1 cMpc, a gas particle mass of 9.2 x 10 4 /i _1 Mq and 
a gravitational softening length of 0.65/i _1 ckpc. The gas is 
reionised instantaneously by a uniform, optically thin ion- 
ising background at z — 9 based on the lHaardt fc Madaul 
(2001) galaxies and quasars emission model. Although we 
subsequently recompute the ionisation state of the gas, this 
means the IGM density field at z < 9 is already Jeans 
smoothed by the increased gas pressure following photo- 
heating (e.g. iPawlik et ahll2009l ). This simulation was pre- 
viously used to model the Lya opacity in the near-zone of 
the high-redshift q uasar ULAS J1120+0641 at 2 = 7.085 
l|Bolton et alj|201ll ). In this work we use an output at z = 7. 

The initial equilibrium ionisation state of the gas is set 
by a spatially uniform ionising background. We do not model 
the patchy ionisation state of the IGM by following the ra- 
diative transfer of ionising photons emitted by individual 
sources. We will, however, estimate the effect of a central 
source on the ionisation state of the IGM later in our anal- 
ysis. A box size of 10h~ cMpc is in any case too small to 
correctly capture the large scale topology of the reionisa- 
tion process. Instead, the high mass resolution of the sim- 
ulation allows us to incorporate the small-scale structures 
responsible for optically thick absorption systems. These ab- 
sorption systems are expected to define the edge of ionised 
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Figure 1. Slices midway (i? x = 5/i _1 cMpc, with width AiJ x = 39ft" 1 ckpc) through the hydrodynamical simulation volume at z = 7 
displaying the spatial distribution of the neutral hydrogen fraction, nm/nn, for four different background photo-ionisation rates. A self- 
shielding prescription (see Eq. [TJ has been used in all four cases. Fully neutral self-shielded regions, which trace the highest density regions 
in the simulation, are shaded white. Anti-clockwise from top left, log^m/s -1 ) = —12.8, —13.6, —14.0 and —13.2. 



bubbles (Crociani ct al. 2011) , regulate the mean free path 
at overlap (IMiralda-Escude et al-lfeood ; iGnedin fc Fanll2006l ; 
iMcQuinn et al.ll201ll) and d elay the latter s t ages of reion- 
isation (|Ciardi et al.1 120061 ; lAlvarez fc Abell 120121 '). These 
systems remain unr esolved in most large-scale reionisation 
simu l ations (but see IChoudhurv et al.l 120091 ; ICrociani et al.l 
120111 : [Al varez fe Abelll2012i for recent "semi-numerical" ap- 



proaches, andlGnedin fc Fanll2006l; iKohler fc Gnedinl [20071 ; 



IMcQuinn et al. 12011 ; Altav et al.ll201ll for smaller simula- 



tions which resolve self-shielded systems) and are therefore 
often not considered when estimating the impact of inter- 
galactic Lya absorption on the visibility of LAEs. 



2.2 Modelling optically thick absorption systems 

In order to include opti cally thick absorbers in our model, 
following ISchavd (|200ll ) we use a simple, physically mo- 
tivated prescription to assign a density threshold above 
which gas re mains self-shielded, which we briefly summarise 
here (see also IMiralda-Escude et al.ll200ol ; iFurlanetto fe Ohl 
2005). Assuming that the typical size of an Hi absorber is 
the Jeans scale and adopting the case-A recombination coef- 



ficient for ionised hydrogen |Abel et al.lll997l ). the overden- 
sity at which an Hi absorber begins to self-shield may be 
approximated as: 



A ss = 36 r 



2/3 , 

12 J 



l + z 



(1) 



where F_i2 = Fhi/10 _12 s _1 is the background photo- 
ionisation rate, T4 = T/10 4 K is the gas temperature, fi is 
the mean molecular weight and f c = n a /nu is the free elec- 
tron fraction with respect to hydrogen. Absorbers at this 
density threshold have a characteristic size 



L ss = 8.7pkpcr:;/ 3 T 4 13 / 30 



V0.61/ V1.08 



1/3 



(2) 



The fiducial values adopted for f e and fj, correspond to a 
plasma of ionised hydrogen and singly ionised helium. Note 
that helium is not expected to be doubly ionised un til z ~ 3 
(e.g. IFurlanetto fc Ohl 120081 ; IMcQuinn et al.l 120091 ). Com- 
parison to radiative transfer calculations indicate this sim- 
ple prescription recovers the density at which intergalac- 
tic gas self-shields to within a factor of two at z = 6 
jMcQuinn et alj|201ll ). 
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Table 1 . The four background photo-ionisation rates used in this 
work. The resulting volume weighted neutral hydrogen fractions 
(which include the contribution from self-shielded gas) are sum- 
marised in the second column. 
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We shall consider four different values for the back- 
ground photo-ionisation rate in this work: log(rm/s _1 ) = 
— 12.8, —13.2, —13.6 and —14.0. These are summarised in 
Table Q] along with the corresponding volume weighted 
Hi fractions. This range is chosen so that the largest value 
is comparable to recent photo-ionis ation rate measurements 
from the Lya forest at z = 6 l|Wvithe fc Bolton] l201ll ; 
ICalverlev et al.|[201ll ). The lowest value produces a volume 
averaged Hi fraction of around ten per cent (with A ss ~ 2), 
and thus approaches the regime where this self-shielding pre- 
scription ceases to provide an adequate description!]] of the 
spatial distribution of H i in the IGM. 

Figure Q] displays the resulting distribution of neu- 
tral hydrogen in a two dimensional slice of the simulation. 
The self-shielding threshold given by Eq. ((TJ is computed 
self-consistently using the temperatures and densities from 
the hydrodynamical simulation. Neutral hydrogen (shaded 
white) traces the densest parts of the cosmic- web, and fills 
a progressively larger fraction of the simulation volume as 
the amplitude of the ionising background is lowered. 



2.3 Extracting line-of-sight quantities 

We use these neutral hydrogen distributions as the starting 
point for constructing mock Lya absorption spectra. Sight- 
lines are extracted parallel to the box boundaries in three 
directions around the 100 most massive dark matter haloes 
in the simulation at z = 7. These halo sight-lines are spliced 
together with sight-lines drawn randomly from the periodic 
volume to form a total of 300 continuous skewers through 
the density, peculiar velocity and temperature field of the 
hydrodynamical simulation. The skewers have a total length 
of 12.5 pMpc, ranging from 8330 kms -1 to -2084kms _1 
around the halo position. The Lya optical depth along each 
of the sight-lines is then co mputed following a standard pro- 
cedure, see e.g. eq. (15) in iBolton fc HaehneltJ i|2007al ). In 
the rest of our analysis, all quantities will be estimated from 
these 300 sight-lines. 

We note, however, that the most massive dark mat- 
ter halo in the simulation volume at z = 7 h a s a m ass 
M dm = 9.0 x 10 9 M Q . In comparison. lOuchi et ail l|2010t ) in- 
fer an average dark matter host halo mass of 10 10 -10 11 M© 
from the clustering of LAEs with L\, ya > 2.5 x 10 42 " 1 



erg s 



1 In other words, when the IGM is no longer in the "post-overlap" 
stage of reionisation (cf. Gncdin 2000) and neutral gas also resides 
in underdense regions far from ionising sources. This self-shielding 
model only applies to what is sometimes referred to as the fi- 
nal "outside-in" stages of reionisation (e.g. Miralda- Escude et al.l 
l200d: IChoudhurv et ' al.ll2009l ; iFinlator et alJl200Sh . 
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Figure 2. Top: comparison of the four photo-ionisation rates 
assumed in this work (filled symbols at z = 7) with obser- 
vational constraints at z < 6 from the Lya forest opacity 
dWvithe fc Boltonl I2011I) and the line-of-sight proximity effect 
llCalverlev et al.l l201l[) . Middle: the filled symbols (offset from 
z = 7 for clarity) display the corresponding volume averaged 
H I fraction. The open symbols show the H I fraction which re- 
sults if self-shielding (see Eq. [TJ is not included. Observational 
constr aints at z < 6.1 are from IGM Lya absorption measure- 
ment s l lFan et al.l2006l;lBolton fc Haehnelt]|2007bl ; [McGreer et~aU 
l201lt) . Bottom: the number density of absorption systems with 
10 17 ' 2 < TVni/cm- 2 < 10 20 per unit redshift. The filled (open) 
symbols at z = 7 again show the model predictions including (ex- 
cluding) gas which is self-shielded from the ionising background. 
The data points with error bars displ ay the number density of ab- 
sorbers per unit redshift measured by[Songaila & Cowie ( 2 010T ) at 
z < 6, rescaled to match the cosmological parameters assumed in 
this work. 



at z — 6.6. Our simulation is therefore too small to correctly 
sample typical LAE host haloes while simultaneously resolv- 
ing the small-scale structure of the IGM, so in what follows 
we will underestimate the impact of the more extended gas 
overdensities and/or stronger inflows expected around larger 
haloes on the visibility of the Lya emitters. However, both of 
these effects will tend to make the impact of Lya absorption 
on the visibility of the Lya emitters stronger, rather than 
weaker, serving to strengthen our main conclusion that mod- 
est Hi fractions may explain the observed LAE fraction at 
z = 7. 



3 PROPERTIES OF THE OPTICALLY THICK 
ABSORBERS 

It is instructive to first make contact with observations by 
comparing our simulation results at z — 7 to constraints 
from the IGM at z < 6. The upper panel in Figure [2] dis- 
plays the four photo-ionisation rates used in Figure [1] The 
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Figure 3. Left: the column density distribution function predicted at z = 7 for the four different background ionisation rates indicated 
on the plot. The vertical dotted lines demarcate the column density thresholds by which we define LLSs (10 17,2 < -/Vjji/cm -2 < 10 19 ), 
super-LLSs (10 19 < 7Vm7 cm ~ 2 < 10 20 3 ) and DLAs (N m > 10 20 3 cm" 2 ). Middle: as for left panel but now displaying log N m f(N m ). 
Right: scatter plots of absorber H I fractions against H I column density for the largest and smallest of the ionisation rates assumed. The 
red triangles show the results obtained when self-shielding is ignored for log^ni/s" 1 ) = —12.8 only. 



filled symbols in the middle panel display the correspond- 
ing volume averaged H i fractions predicted by these models. 
The lowest Hi fraction is consistent with only weak evolu- 
tion from the lower observational limit at z = 6, whereas 
smaller photo-ionisation rates result in a stronger increase 
in the Hi fraction toward z — 7. Note these Hi fractions 
are higher than the values expected if self-shielding were ig- 
nored - this is demonstrated by the open symbols at z ~ 7, 
which display the H i fraction computed from our simulated 
sight-lines without the self-shielding prescription given by 
Eq. {1} (i.e. assuming a uniform, optically thin ionising back- 
ground) . 

The bottom panel of Figure [2] displays the resulting 
number density of absorption systems with column^ den- 
sities 10 17 - 2 < 7V H i/cm~ 2 < 10 20 , alo ng with the recent 
observ ational measurements presented bv lSongaila fc Cowid 
(2010) at z < 6. The filled (open) symbols again display 
the results including (excluding) self-shielding. Lowering the 
amplitude of the Hi ionising background significantly in- 
creases the incidence of these absorption systems. This in- 
crease is more pronounced in the models incorporating self- 
shielding, with a factor of ~ 25 increase in the number 
density of absorbers per unit redshift from log(FHi/s _1 ) = 
-12.8 to -14.0. 

The absorber properties are displayed in more detail 
in Figure [3] where the left panel displays the H i column 
density distribution functions (CDDFs) obtained from the 
simulations. The CDDF follows the conventional definition: 



d 2 N 



HI, Z 



H(z) 



dzdNm H (l + z) 



(3) 



Including self-shielding increases the number of optically 
thick absorbers, flattening the otherwise power-law column 



2 Column densities are estimated by integrating the density field 
over a scale of 20 pkpc; this choice is motivated by Eq. (0, which 
yields L ss ~ 16-40 pkpc for the range of photo-ionisation rates 
considered here. 



density distribution at 10 17 ' 2 < AT H i/cm~ 2 < 10 19 ' 2 . A 
similar flattening of the CDDF is inferr ed from observa- 
tions at z ~ 3.7 |Prochaska et al]|2010t ). and it has also 
been noted in many other theoretical models jKatz et al.l 
199fj; IZheng fc Miralda-Escudeil2002l; iNagamine et al.ll201of 



McQuinn et al.ll201ll ; lAltav et all l201ll ). The distributions 
converge at Am ~ 10 20 cm~ 2 , indicating the highest col- 
umn density systems (damped Lya absorbers with Am > 
10 20 ' 3 cm -2 ) remain largely insensitive to changes in the 
background photo-ionisation rate by virtue of their high gas 
density, at least for Thi < 10~ 12,8 s -1 . The middle panel 
in Figure [3] displays Ahi/(Ahi), which indicates which op- 
tically thick systems comprise the bulk of the Hi opacity. 
When self-shielding is included, it is absorption systems with 



At 



10 1 



-10 1 



' cm which dominate. 



Finally, the ionisation state of the optically thick ab- 
sorbers ar e displayed in the righ t panel of Figure (see also 
figure 5 of lMcQuinn et aI1l201ll ) . The blue diamonds and or- 
ange squares show the average Hi fraction against column 
density for log^m/s" 1 ) = —12.8 and —14.0, respectively. 
As the photo-ionisation rate is lowered, the density thresh- 
old for self-shielding is reduced (since A ss oc r_-f 2 ) and the 
size and incidence rate of optically thick absorption systems 
(i.e. Am > o^ 1 = 10 17 ' 2 cm -2 , where ohi is the Hi photo- 
ionisation cross-section at the Lyman limit) increases. The 
self-shielded regions first extend into the outer parts of galac- 
tic haloes and then into the filaments defining the cosmic 
web (e.g. Figure [1}. For comparison, the red triangles dis- 
play the results obtained for log^m/s -1 ) = —12.8 when 
ignoring self-shielding. Absorbers with the same A^hi have 
larger Hi fractions when self-shielding is included, imply- 
ing these must correspond to regions of lower gas density 
relative to the optically thin case. As we shall now demon- 
strate, it is this rapid increase in the incidence of optically 
thick absorption systems with increasing neutral fraction at 
the tail-end of reionisation which is primarily responsible for 
strengthening the red Lya damping wing in our simulations. 
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Figure 4. The median Lya transmission redward of the rest frame Lya transition (at Av = Okms -1 ) obtained from an ensemble of 300 
simulated sight-lines. Each panel displays a different intergalactic Ly-cr absorption model. The four continuous curves in each panel display 
the transmission profile for different background ionisation rates as indicated. The binned grey line displays the observed transmitted 
fraction around the z = 7.085 quasar ULAS J1120+0641 with an estimated 7V ion = 1.3 X 10 57 s~ 1 . Left: model which includes self-shielded 
(SS) gas, but ignores the influence of a central ionising source. Centre: as for left panel, but now each sight-line includes the photo-ionisation 
rate from an LAE with LLya = 5 X 10 42 ergs -1 which emits ionising photons at a rate JVi on = 1.1 X 10 53 s -1 . Right: a model which neglects 
optically thick systems, and instead models the ionisation state of the IGM as an ionised bubble with radius Run embedded in an otherwise 
fully neutral IGM. The Hll bubble sizes are as indicated on the plot. Note this case is shown to emulate results from existing models in 
the literature only. 



4 IMPLICATIONS FOR LAE VISIBILITY 

4.1 Photo-ionisation by a central source 

We now turn to the main result in this paper; estimating 
the impact of the increasing incidence of optically thick ab- 
sorbers on the visibility of LAEs. Before proceeding further, 
however, we also allow for the possible effect of ionising radi- 
ation from a central source in our models (i.e. the Lya emit- 
ter under observation). Ionising radiation from the LAE it- 
self can reduce the incidence of proximate self-shielded gas 
and therefore weaken the Lya damping wing. 

We model the possible effect of the LAE as follows. 
Firstly, we compute the emission rate of ionising photons 
into the IGM, Ni on [s _1 ], by an LAE with observed Lya lu- 
minosity LLya assuming case-B recombination: 



AT — ^ i'Lya 



fm 



2 hpV a (1 — fesc)fa 



(4) 



Here the energy of a Lya photon is h-pv a — 10.2 eV, 
and /esc and f a are the escape fractions of Lyman con- 
tinuu m and Lya ph o tons into the IGM. In this work we 
follow iDiikstra et al. | (|2007h and assume f a = 1, such that 
Lhya = / a Ljf'* yields the Lya luminosity uncorrected for 
dust. We adopt Z/Lya = 5 x 10 42 ergs _1 based on the typ- 
ical Lya luminos ities measured in rec ent observations of 
z ~ 7 LAEs (e.g. IPentericci et al.ll201ll ). We assume a Ly- 
man continuum escape fraction of / csc = 0.2, which yields 



A^ion = 1.1 x 10° 



The Lyman continuum escape frac- 



tion is highly uncertain. A value > 20 per cent is required 
at z ~ 6 if the observed Lyman break galaxy popula- 
tion is to repr oduce ionising background constraints from 
the Lya forest |Bolton fc Haehnelt 2007b; S hull et alJlioil; 



iFinkelstein et al.ll2012T ). but as we will discuss in more de- 
tail in Section 14.41 the relation between Lya emission and 
the escape of ionising radiation may not be straightforward. 

In the models with a central ionising source we then 
assume the total photo-ionisation rate to be the sum of the 
LAE contribution and the metagalactic background rate: 



CHI 



JVion g 

AttR 2 /3 + 3 



(5) 



where we have assumed a power law spectral energy distri- 
bution for the io nising radiation from the LAE, L v oc v~^, 
with p = 3 (cf. iLeitherer et~al]|l999h . Note that we have 
ignored (i) the finite lifetime of the LAE (ii) the inter- 
vening Hi opacity between self-shielded systems and (iii) 
assumed the LAE luminosity is constant. All three as- 
sumptions will overestimate the impact of ionising radia- 
tion from the LAE on the surrounding neutral hydrogen. In 
what follows we therefore consider this as a limiting case 
and bracket other possibilities by additionally computing 
the Lya damping wing when the LAE is ignored. We do 
not explicitly account for the possible clustering of ionising 
sources around the LAE in addition to the range of back- 
ground photo-ionisation rates we already consider. As we 
shall demonstrate in the next section, larger background 
photo-ionisation rates weaken the red Lya damping wing 
by raising the density threshold at which gas self-shields. 



4.2 The red Lya damping wing 

Figure 2] displays the median IGM Lya transmission pro- 
file obtained from the 300 simulated sight-lines - each panel 
corresponds to a different IGM absorption model. Note that 
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when computing the transmission profiles we have ignored 
the contribution of all gas within 20 pkpc of the centre of the 
LAE host halo to both the Lyman continuum and Lya opac- 
ity. We do not model the complex radiative transfer within 
the host galaxy's interstellar medium, and account for the 
absorption by this material with the parameters / esc and f a 
instead. 

The left panel shows the transmission profile predicted 
by our simulations without photo-ionisation by the LAE. A 
strong damping wing profile is produced, even for photo- 
ionisation rates that correspond to Hi fractions as low as 
~ 3 per cent (blue dot-dashed curve). This is due to the 
presence of high column density LLS or damped Lya ab- 
sorption systems within ~ 1 pMpc of the LAE host halo 
(at Av — Okms -1 ). For comparison, the grey binned line in 
all panels displays the tran smission profile of the z = 7.085 
quasar ULAS J1120+0641 (|Mortlock et al.ll201ll ). 

The central panel instead shows the damping wings for 
our models including ionising emission from an LAE. The 
source dominates the total photo-ionisation rate (i.e. pro- 
duces more than 50 per cent of the total) within ~ 0.5 pMpc. 
The incidence of proximate, self-shielded systems within this 
radius is therefore significantly reduced, and the Lya damp- 
ing wing is consequently weakened. On the other hand, re- 
call that the typical halo mass in our simulations is rather 
small (see Section |2.3[) . We therefore likely underestimate 
the abundance of proximate optically thick systems and 
this effect may be somewhat less important. The conver- 
gence of the transmission profiles at F — (close to Av = 
— 50 kins -1 ) is due to peculiar velocities arising from grav- 
itational infall around the haloes. We have verified that ig- 
noring peculiar motions produces profiles which converge 
at Av — Okms -1 instead. Note that in the left panel of 
Figure 21 where emission from a central LAE is excluded, 
this convergence is less prominent. This is because absorp- 
tion line cores from proximate absorbers in addition to the 
Lya damping wing can extend redward of systemic Lya; 
this occurs more frequently as the photo-ionisation rate is 
lowered and the size of the absorbers increases. 

At this point, it is important to emphasise that the 
transmission profiles also vary significantly between individ- 
ual sight-lines; a strong damping wing will occur whenever 
there are super-LLS or damped Lya absorbers lying close 
to the host halo. This can occur even when the IGM is 
highly ionised (i.e. for log(rHi/s -1 ) = —12.8) assuming a 
sufficiently overdense region lies in close proximity to the 
LAE. This variation can be seen clearly in Figure [5l which 
shows the 68 and 95 per cent bounds for the scatter around 
the median Lya absorption profile for log(FHi/s _1 ) = —12.8 
(upper panel) and —14.0 (lower panel). Some Lya emitters 
will therefore be completely obscured by the damping wing 
while others will remain fully visible. The stochastic nature 
of the damping wing absorption from these high column den- 
sity systems should therefore also impact on the clustering 
of LAEs (e.g. iFurlanetto et alJliooil ; iMcQuinn et alj|2007al ) 
for even relatively low H i fractions of 3-10 per cent, mod- 
ifying the enhanced clustering signal expected as the IGM 
Hi fraction rises. This will be an interesting possibility to 
explore in future work. 

In summary, it is clear optically thick absorbers can 
produce a rapid decline in the transmission of Lya emis- 
sion from high-redshift galaxies at the tail-end of reioni- 
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Figure 5. The scatter in Lya transmission along individual sight- 
lines. Each panel assumes the background photo-ionisation rates 
indicated on the diagram and include a central ionising source 
with N ion = 1.1 X 10 63 s — The short-dashed (red) and long- 
dashed (blue) curves bound 68 and 95 per cent of the Lya trans- 
mission around the median (black curves). Note the median pro- 
files are identical to those displayed in the central panel of Fig- 
ure [4] 



sation. This is quite different to existing predictions us- 
ing large-scale numerical and semi-analytical simulations of 
reionisation. In the late stages of reionisation the presence 
of self-shielded regions where h ydrogen can recombine and 
stay neutral is ver y important JlVIiralda-Escude et al.l 2000l; 
Ciardi et al l l200d iGnedin fc Fan! 120061; IChoudhurv et al l 



20091 ; ICrociani et alj l201ll ; lAlvarez fc Abell |2012( ). Large 
scale reionisation simulations do not yet h ave the dynamic 
range to resolve these systems (see e.g. iTrac fc Gnedinl 
1201 ll ). Any neutral gas is thus typically far away from the 
central regions of Hn regions, and the impact of nearby op- 
tically thick absorbers on the Lya damping wing is ignored. 

We illustrate this in the right-hand panel of Figure 2J 
where we have simulated the Lya damping wing expected 
around a LAE embedded in an H n bubble with radius Rhii. 
Inside the bubble (R < -Rhii), the IGM is highly ionised with 
a uniform H 1 fraction of 10 -4 , while outside the IGM is fully 
neutral. This simple model serves to illustrate the impor- 
tance of optically thick absorbers for the red damping wing. 
Bubble sizes of Rhii < 30 cMpc are required to produce a 
damping wing similar in strength to the absorption seen 
when the IGM is only ~ 3 per cent neutral by volume (blue 
dot-dashed curve in the central panel of Figure |4} when in- 
cluding optically thick absorbers. In comparison, large scale 
reionisation simulations typically predict Hn bubble sizes 
of ~ 30-40 cM pc for Hi fraction s as large as 30 per cent 
(e.g. figure 2 in lZahn et al.1120111 ). This elucidates why ex- 
isting models indicate the LAE/LBG fraction at z ~ 7 may 
imply large Hi fractions of 40-90 per cent. At fixed neutral 
fraction, large-scale reionisation simulations predict a damp- 
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ing wing which is significantly weaker compared to the case 
where absorption by optically thick absorbers is included. 
The damping wing will therefore only reduce LAE visibility 
at times when individual H n bubbles have not yet grown to 
the typical size of > 50 cMpc reached at overlap in these 
simulations. 



4.3 The REW distribution 

We now compare our models directly t o recent mea- 
surements of the LBG/LAE fraction 



(Font ana ct al 

2010| : lyanzella et al.ll2oTTI:|Pentericci et al-lfeoilUOno et al 



2012 



Schcnkcr ct al. 



I l2012h 



following the approach of 



Diikstra et alJ ( 201 lh . which we recapitulate here. Our first 



requirement is to compute the fraction of the Lya emission 
line which is transmitted, SF ', where 



,9 = 



(6) 



Here J{v) gives the shape of the Lya line profile, e _TLyQ ^- ) 
is the IGM transmission (e.g. Figure [J) and i-w n and v max 
are chosen to correspond to the full extent of our simulated 
sight-lines. To compute ST for each of our simulated sight- 
lines we must assume an intrinsic Lya profile. We adopt a 
Gaussian line profile with a standard deviation of 88 kms -1 , 
corresponding to the circular veloci ty of a dark matter hal o 
of mass 10 10 M Q at z = 7 (eq. 25 in lBarkana fc LoebllioOll ). 
We further assume that the line centre is shifted redward 
of Lya by 400 kms -1 ; this offset is typi cal of Lya emissio n 
lines associated with LBGs at z ~ 2 — 3 (St eidel et al.ll2010f ). 
and is chosen to mimic the effect of resonant scattering by 
the (outflowing) gas on the Lya line visibility (but see our 
discussion in the next section). 

We next construct a probability distribution dP(< 
log S/"} I d log 3? for each IGM absorption model using our 
simulations. For Lya redshifted relative to systemic we as- 
sume the observed rest frame equivalent width (REW) at 
2 = 6, when the hydrogen in Universe is already highly 
ionised, to be equal to the emitted REW (i.e. 3T z =q = 1). If 
we further assume that there is no intrinsic evolution in the 
LAE population from z — 6 to z — 7, the cumulative prob- 
ability distribution (CPDF) that a galaxy has an observed 
REW and an IGM transmission ST at z — 7 is 



P(> REW) 



d „ dP{< logf) P, =6 (> REW/.^) 



dlog^ 



^lnlO 



Here REWi„t r = REW/, 9, P z = 6 (> REW) = 
exp(- REW/REW c ) and REW C = 50 A (see iDiikstra et all 
l201ll for details). 

Figure \E\ displays the resulting CPDFs for the IGM 
models with self-shielding only (left panel) and including 
ionisation by a central LAE (central panel). The CPDFs 
are compa r ed to a recent data compilation presented by 
Ono et all ||2012[) . as well as an upper limit at 75 A from 



Stark etaiM . As expected, our models predict a rapid 



evolution in the REW distribution from z = 6, even for mod- 
est neutral fractions of 3 per cent (blue dot-dashed curve). 
In the right panel, the CPDF for the log(r H i/s -1 ) = -14.0 
model in the central panel (orange long-dashed curve with 
(/hi)v ~ 0.1) is shown assuming four different offsets for 
the Lya emission line relative to systemic. Shifting the line 



blueward (redward) of the fiducial —400 kms -1 decreases 
(increases) producing a stronger (weaker) evolution in 
the REW distribution. 



4.4 Uncertainties 

From our discussion thus far, it should be clear that mod- 
elling the IGM transmission for LAEs at the tail-end of 
reionisation is uncertain mainly because of the difficulties 
with self-consistently modelling self-shielded regions. How- 
ever, uncertainties associated with the escape of ionising ra- 
diation from the (possibly clustered) Lya emitters and the 
redshift relative to systemic of the Lya emission also play a 
role. 

When including an LAE in our models we have as- 
sumed / csc = 0.2. It is, however, not obvious that high red- 
shift galaxies emit ionising and Lya radiation at the same 
time and/or in the same direction. It is certainly plausible 
that no ion ising radiation esca pes in the direction of the ob- 
server (e.g. iGnedin et al]|2008l ), in which case even our IGM 
model excluding a central source may actually be an overes- 
timate of the IGM transmission 3?\ recall we have ignored 
neutral gas within 20pkpc of the halo centre and instead 
assumed the appropriate escape fractions for ionising and 
Lya photons. Note further that the escape of ionising ra- 
diation and the emission redshift relative to systemic may 
also be unfavourably coupled. A large redshift for the emis- 
sion line may need a large column density of neutral gas 
within the host halo of the LAE. This is required in order to 
resonantly scatter the Lya photons local ly far into the red 
wing of the line (e.g. iBarnes et alJl201lT l. The redshifts of 
400 km s -1 or more for the observed Lya emission in z = 3 
LBGs (|Steidel et alJ l2010h may therefore be incompatible 
with a simultaneous escape of ionising radiation along the 
line-of-sight to the observer. 

The detailed intrinsic line Lya shape is also rather un- 
certain - in this work we have assumed a simple Gaussian 
line profile with standard deviation 88 kms -1 . In practice, 
however, radiative transfer models predict a characteristic 
double peaked emission line (e.g. IVerhamme et all 120061 ) 
with intergalactic Lya absorption removing the blue peak 
at high redshift. In their recent analysis of the REW dis- 
tribution predicted by large-scale reionisation simulations, 
Ijensen et al.l (|2012T ) adopt a doubled peaked Gaussian mi- 
nus a Gaussian line based on the r adiati ve transfer simu- 
lations presented by lLaursen et afl (|201ll ). They note that 
using this line profile weakens the impact of the H 1 fraction 
on the REW distribution compared to a simple Gaussian 
centred at the systemic redshift. Note, however, since we 
adopt a profile which is shifted redward of systemic, our as- 
sumption that the entire line is visible in the absence of a 
red Lya damping wing is reasonable. Unless the line width 
becomes comparable to the line offset relative to systemic, 
this will not have a significant effect on the resulting REW 
distribution. 

Lastly, the small simulation box size of 10/i -1 cMpc and 
the omission of self-consistent modelling for the ionising 
sources remains the most significant concern in this anal- 
ysis. We do not correctly model large-scale variations in the 
IGM ionisation state prior to overlap. On the other hand, we 
argue that the proper modelling of dense, self-shielded re- 
gions are an equally important part of this problem. A fully 
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Figure 6. Predicted cumulative probability distributions (CPDFs) for the rest frame Lya equival ent width (REW) at 2 = 7. The dotted 
curves display the intrinsic 2 = 6 CPDF assuming the IGM is fully transparent to Lya photons jDiikstra et al.ll201lT ). The solid, short- 
dashed, dot-dashed and long-dashed curves are constructed from our Lya transmission models at 2 = 7 excluding (left) and including 
(centre) ionising radiation emitted along the line of sight from a source with Z/Lya = 5 X 10 42 ergs -1 and N[ on = 1.1 X 10 53 s — 1 (see text for 
details). These models adopt an intrinsic REW distribution identical to that at 2 = 6 and differ due to evolution in the Lya transmission 
only. A Gaussian Lya profile with a standard deviation of 88kms -1 and an offset of — 400kms~ 1 redward of rest frame Lya has been 
assumed in both cases. The right hand panel instead displays the log(rHi/s — 1 ) = —14.0 case from the central pa nel (with (/hi)v ~ 0.1) 
assuming four different Lya emission lin e offsets. The models a re compared to t he composite data set c o m| 3iled bv lOno et all i2012l) (fille d 
circles) which utilises observations from iFontana et al.l lj2010l ) ; |Vanzelja^tal 1 ll201ll ); |p"entericci et"aH l|201ll ) and ISchenker et al.l l|2012l ). 
The upper limit at 75 A is based on the earlier non-detection bv lStark et al.1 J20ld) at 2 ~ 7. 



self-consistent model will require combining both of these 
aspects, either by resolving all relevant physical scales or 
adopting sub-grid models which are informed by high resolu- 
tion simulations of the IGM density field. Furthermore, any 
modelling should be calibrated to reproduce the observed 
transmission redwards of Lya in bright high-redshift quasars 
which appears to evolve rapidly between z = 6-7; the simu- 
lations presented here have been used to estimate a neutral 
fraction of > 10 p er cent around the z = 7.085 quasar ULA S 
J1120+0641 (e.g. iMortlock et alJlioill ; iBolton et al.ll201ll ). 
Note, however, this is based on a small number of bright 
quasars at z ~ 6 and only one at z > 7. Future observations 
will help determine whether ULAS J1120+0641 is typical. 



| 2010l ; iPentericci et all l201ll; lHaves et al.1 l201ll; lOno et al.1 
hoid i lSchenker et al.ll2012l ; ICurtis-Lake et alj|2012ft . We find 
that these observations do not require a large neutral frac- 
tion (~ 40-90 per cent) in the intervening IGM at z = 7 
as previously suggested. Instead, if the rapid decline in the 
LAE/LBG fraction is further corroborated, it may instead 
be indicative of the rapid decrease of the mean free path 
of ionising photons expected at the tail-end of reionisation. 
Furthermore, as we find the Lya emission from high-redshift 
galaxies will be suppressed for volume-averaged neutral frac- 
tions of only 3-10 per cent, the patchiness of reionisation 
may produce a more modest impact on the clustering prop- 
erties of observable high-redshift LAEs than previously pre- 
dicted. 



5 CONCLUSIONS 

We have used a hydrodynamical simulation to model the 
Lya opacity of the intervening IGM during the final stages 
of reionisation. As the photo-ionisation rate drops the opac- 
ity redward of rest-frame Lya is expected to rise rapidly due 
to the increasing incidence of optically thick absorption sys- 
tems. Our results indicate that the bulk of the H 1 opacity 
will arise from optically thick systems with column densities 
N m ~ 10 18 ' 5 -10 19 ' 5 cm' 2 . When including these absorption 
systems in the simulations, only a moderate rise in the vol- 
ume averaged neutral fraction (to (/hi)v ~ 0.03-0.1) is re- 
quired to significantly reduce the transmission redward of 
rest-frame Lya. 

This result has an important implication for the inter- 
pretation of the recently observed decline in th e Lya emis- 
sion from high-redshift galaxies at z ~ 6-7 (e.g. IStark et al.l 



Our findings may be particularly relevant for future sur- 
veys which plan to use the transmission of Lya emission and 
the clustering properties of LAEs to probe deep into the 
epoch of reionisation. Detecting these LAEs at z > 7 may 
well become difficult, even when the IGM is only ten per cent 
neutral. It may also mean that spectroscopic confirmation of 
high-redshift candidate galaxies identified with the drop-out 
technique at z > 7 may become problematic, at least un- 
til emission lines other than Lya can be observed with the 
James Webb Space Telescope. Finally, extending theoreti- 
cal predictions beyond the approximate approach adopted 
in this work is vital. This will ultimately require incorporat- 
ing small-scale absorption systems into simulations which 
also model the patchy nature of reionisation on large scales. 
This remains a considerable but fundamental computational 
challenge for the current generation of reionisation simula- 
tions. 
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